// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                  avtAggregateChordLengthDistributionQuery.C               //
// ************************************************************************* //

#include <avtAggregateChordLengthDistributionQuery.h>

#include <vtkCellData.h>
#include <vtkIdList.h>
#include <vtkPolyData.h>

#include <avtLineScanFilter.h>
#include <avtParallel.h>
#include <avtOriginatingSource.h>

#include <DebugStream.h>

#include <stdio.h>
#include <math.h>

#include <vector>

// ****************************************************************************
//  Method: avtAggregateChordLengthDistributionQuery constructor
//
//  Purpose:
//      Defines the constructor.  Note: this should not be inlined in the header
//      because it causes problems for certain compilers.
//
//  Programmer: Hank Childs
//  Creation:   August 31, 2006
//
// ****************************************************************************

avtAggregateChordLengthDistributionQuery::avtAggregateChordLengthDistributionQuery()
{
    numChords = new int[numBins];
}


// ****************************************************************************
//  Method: avtAggregateChordLengthDistributionQuery destructor
//
//  Purpose:
//      Defines the destructor.  Note: this should not be inlined in the header
//      because it causes problems for certain compilers.
//
//  Programmer: Hank Childs
//  Creation:   August 31, 2006
//
// ****************************************************************************

avtAggregateChordLengthDistributionQuery::~avtAggregateChordLengthDistributionQuery()
{
    delete [] numChords;;
}


// ****************************************************************************
//  Method: avtAggregateChordLengthDistributionQuery::PreExecute
//
//  Purpose:
//      Does some initialization work before the query executes.
//
//  Programmer: Hank Childs
//  Creation:   August 31, 2006
//
// ****************************************************************************

void
avtAggregateChordLengthDistributionQuery::PreExecute(void)
{
    avtLineScanQuery::PreExecute();

    delete [] numChords;
    numChords = new int[numBins];
    for (int i = 0 ; i < numBins ; i++)
        numChords[i] = 0;
}


// ****************************************************************************
//  Method: avtAggregateChordLengthDistributionQuery::PostExecute
//
//  Purpose:
//      Outputs the chord counts.
//
//  Programmer: Hank Childs
//  Creation:   August 31, 2006
//
//  Modifications:
//    Eric Brugger, Thu Mar 25 09:42:25 PDT 2010
//    I changed the name of the curve in the ultra file to avoid using
//    special characters.
//
//    Kathleen Biagas, Tue Feb 25 08:58:52 PST 2014
//    Add XML results, and ResultValues, allowing them to be set even if
//    output file could not be opened.
//
// ****************************************************************************

void
avtAggregateChordLengthDistributionQuery::PostExecute(void)
{
    int   i;

    int times = 0;
    char name[1024];
    sprintf(name, "cld_a%d.ult", times++);

    if (PAR_Rank() == 0)
    {
        bool lookingForUnused = true;
        while (lookingForUnused)
        {
            ifstream ifile(name);
            if (ifile.fail())
                lookingForUnused = false;
            else
                sprintf(name, "cld_a%d.ult", times++);
        }
    }

    char msg[1024];
    sprintf(msg, "The aggregate chord length distribution has been outputted "
            "as an Ultra file (%s), which can then be imported into VisIt.", 
            name);
    SetResultMessage(msg);
    SetResultValue(0.);

    int *nc2 = new int[numBins];
    SumIntArrayAcrossAllProcessors(numChords, nc2, numBins);
    delete [] numChords;
    numChords = nc2;

    if (PAR_Rank() == 0)
    {
        double binWidth = (maxLength-minLength) / numBins;
        double totalArea = 0.;
        for (i = 0 ; i < numBins ; i++)
            totalArea += binWidth*numChords[i];
        if (totalArea == 0.)
        {
            sprintf(msg, "The chord length distribution could not be "
                    "calculated because none of the lines intersected the data"
                    " set.  If you have used a fairly large number of lines, "
                    "then this may be indicative of an error state.");
            SetResultMessage(msg);
            return;
        }

        ofstream ofile(name);
        if (ofile.fail())
        {
            sprintf(msg, "Unable to write out file containing distribution.");
            SetResultMessage(msg);
        }
        if (!ofile.fail())
            ofile << "# Chord length distribution - aggregate" << endl;

        MapNode result_node;
        doubleVector curve;

        for (int i = 0 ; i < numBins ; i++)
        {
            //double x = minLength + (i+0.5)*binWidth;
            double x1 = minLength + (i)*binWidth;
            double x2 = minLength + (i+1)*binWidth;
            double y = numChords[i] / totalArea; // Make it be 
                            // a distribution ... the area under the curve: 1
            curve.push_back(x1);
            curve.push_back(y);
            curve.push_back(x2);
            curve.push_back(y);
            if (!ofile.fail())
            {
                ofile << x1 << " " << y << endl;
                ofile << x2 << " " << y << endl;
            }
        }
        result_node["chord_length_distribution_aggregate"] = curve;
        SetXmlResult(result_node.ToXML());
        SetResultValues(curve);
    }
}


// ****************************************************************************
//  Method: avtAggregateChordLengthDistributionQuery::ExecuteLineScan
//
//  Purpose:
//      Examines the input data.  Note that the line scan filter will organize
//      the data so that all of the contributions from a given line will be in
//      the same vtkPolyData input.
//
//  Programmer: Hank Childs
//  Creation:   August 31, 2006
//
// ****************************************************************************

void
avtAggregateChordLengthDistributionQuery::ExecuteLineScan(vtkPolyData *pd)
{
    vtkIntArray *lineids = (vtkIntArray *) 
                                  pd->GetCellData()->GetArray("avtLineID");
    if (lineids == NULL)
        EXCEPTION0(ImproperUseException);
        
    int npts = pd->GetNumberOfPoints();
    std::vector<bool> usedPoint(npts, false);
    
    pd->BuildLinks();
    pd->BuildCells();

    // When we determine which line segments are on one side of another,
    // we need to examine the other segments.  But we only want to consider
    // segments that are from the same line (i.e. have the same lineid).
    // So we are using this data structure to keep track of which segments
    // come from which line.  (Instead of having to iterate through all the
    // other segments each time we examine a segment.)
    //
    // So why hash?
    // The lines are spread out over many processors.  So proc X may have
    // lines 100,000 - 100,500.  So we either have to have vectors that are
    // way too big (100,500 entries with the first 100,000 empty),
    // or we have to do some clever indexing.  So we are doing
    // some clever indexing.  In this case, hashing.
    int hashSize = 10000;
    std::vector< std::vector<int> >    hashed_lineid_lookup(hashSize);
    std::vector< std::vector<double> > hashed_segment_length(hashSize);

    for (int i = 0 ; i < npts ; i++)
    {
        if (usedPoint[i])
            continue;
        int seg1, seg2;
        int numMatches = GetCellsForPoint(i, pd, lineids, -1, seg1, seg2);
        if (numMatches == 0)
            continue;
        if (numMatches > 2)
        {
            // We found an error condition.  Give up on this line.  This
            // happens infrequently enough that it should not affect our
            // statistics.
            continue;
        }
        int oneSide = i;
        int otherSide = i;
        int lineid = lineids->GetValue(seg1);
        if (numMatches == 1)
        {
            oneSide   = i;
            otherSide = WalkChain(pd, i, seg1, usedPoint, lineids, lineid);
        }
        else if (numMatches == 2)
        {
            oneSide   = WalkChain(pd, i, seg1, usedPoint, lineids, lineid);
            otherSide = WalkChain(pd, i, seg2, usedPoint, lineids, lineid);
        }
        if (oneSide == -1 || otherSide == -1)
        {
            // We found an error condition.  Give up on this line.  This
            // happens infrequently enough that it should not affect our
            // statistics.
            continue;
        }
        double pt1[3];
        double pt2[3];
        pd->GetPoint(oneSide, pt1);
        pd->GetPoint(otherSide, pt2);
        double dist = sqrt((pt2[0]-pt1[0])*(pt2[0]-pt1[0]) + 
                           (pt2[1]-pt1[1])*(pt2[1]-pt1[1]) + 
                           (pt2[2]-pt1[2])*(pt2[2]-pt1[2]));
 
        int hashid = lineid % hashSize;
        hashed_lineid_lookup[hashid].push_back(lineid);
        hashed_segment_length[hashid].push_back(dist);
    }

    for (int i = 0 ; i < hashSize ; i++)
    {
        std::vector<int> already_considered;
        for (size_t j = 0 ; j < hashed_lineid_lookup[i].size() ; j++)
        {
             bool alreadyDoneLineId = false;

             for (size_t k = 0 ; k < already_considered.size() ; k++)
                 if (hashed_lineid_lookup[i][j] == already_considered[k])
                     alreadyDoneLineId = true;
             if (alreadyDoneLineId)
                 continue;

             int lineid = hashed_lineid_lookup[i][j];
             already_considered.push_back(lineid);
             double length = 0.;
             for (size_t k = j ; k < hashed_lineid_lookup[i].size() ; k++)
             {
                 if (hashed_lineid_lookup[i][k] == lineid)
                     length += hashed_segment_length[i][k];
             }

             int bin = (int)((length-minLength)/(maxLength-minLength) * numBins);
             if (bin < 0)
                 bin = 0;
             if (bin >= numBins)
                 bin = numBins-1;
             numChords[bin]++;
        }
    }
}


